Genotype-by-environment interaction and genetic dissection of heartwood color in Cryptomeria japonica based on multiple common gardens and quantitative trait loci mapping

The heartwood color of a major plantation tree Cryptomeria japonica shows high variability among clones and cultivars, and brighter heartwood has higher value in the usage of non-laminated wood such as in traditional construction, which makes heartwood color an important trait in breeding of this species. However, the genetic basis of the interactions between genetics and the environment on heartwood color has been understudied while these are necessary for effective breeding programs in multiple environmental condition. The objectives of the present study were to evaluate the effects of genetics and environments on heartwood color and how they interact in contrasting environments, and to identify genomic regions controlling heartwood color in C. japonica across multiple environments. Heartwood color in terms of L*a*b* color space and spectral reflectance was measured in common gardens established in three contrasting sites. Quantitative trait loci (QTL) that affect heartwood color were identified using previously constructed highly saturated linkage maps. Results found that heartwood color was largely genetically controlled, and genotype-by-environment interaction explained one-third of the total genetic variance of heartwood color. The effect of the environment was small compared to the effect of genetics, whereas environmental effects largely varied among heartwood color traits. QTL analysis identified a large number of QTLs with small to moderate effects (phenotypic variation explained of 6.6% on average). Some of these QTLs were stably expressed in multiple environments or had pleiotropic effects on heartwood color and moisture content. These results indicated that genetic variation in phenotypic plasticity plays an important role in regulating heartwood color and that the identified QTLs would maximize the breeding efficiency of heartwood color in C. japonica in heterogeneous environments.


Introduction
In forest trees, a number of studies on the inheritance and dissection of quantitative traits have focused on wood properties and growth (reviewed in [1]). In wood property traits, strength and quality-related traits, such as Young's modulus and density, have been studied due to their importance in wood characteristics [2][3][4]. However, breeding of quantitative traits in tree species remains challenging because of long generation time of forest trees. Previous studies showed that selection for traits of interests based on marker-based genotypes which is known as marker-assisted selection (MAS) and selection of genotypes based on genomic breeding values derived from genome-wide DNA markers (genomic selection) substantially enhance breeding efficiency, and breeding approaches have attracted interest of plant breeders [5][6][7].
Quantitative traits, such as growth and wood properties, are characterized by genetic and environmental factors and interactions between genetics and environments [8][9][10]. Understanding genetic and environmental effect changes and genotype-by-environment interaction (genetic variation of phenotypic plasticity) in different environments is of primary importance in forest tree breeding and conservation under environmental changes, such as rapid climate changes and seed/seedling transfer via assisted migration [11]. Common garden experiments in multiple environments have provided important insights into the relative contribution of genetic and environmental factors on phenotypic traits, and genes or genomic regions that affect phenotypic traits have been identified using quantitative trait locus (QTL) and association mapping in forest trees (reviewed in [12]). Accumulation of the knowledge of QTLs would contribute to breeding of quantitative traits, however, the effect of QTLs could also change depending on the environmental conditions [13]. Thus, understanding the magnitude of these interactions across multiple environments is necessary for understanding the genetic architecture and improving the breeding efficiency of quantitative traits.
Wood color is one of the traits that most differentiate wood properties and is considered to enhance the value of the wood. In this context, wood color variations have been studied in various tree species [14][15][16]. However, studies on wood color trait dissection have remained scarce compared to other wood properties, such as wood strength. In some forest tree species, wood color is highly correlated with other wood properties, such as wood density and moisture content [14,[17][18][19]. For example, Rink and Phelps [14] reported a negative correlation between wood density and color (lightness) in Juglans nigra. Moya et al. [18] reported that wood density was negatively correlated with moisture content and b � (blue/yellow) color in teak (Tectona grandis). The esthetic value and large variability of heartwood color, and the strong correlations of heartwood color with other wood property traits described above imply that heartwood color is suitable for tree breeding and it may be possible to improve breeding efficiency of heartwood color through molecular breeding approaches.
Cryptomeria japonica (D. Don) is an allogamous, wind-pollinated, and evergreen conifer species [20]. This species is one of the most commercially important tree species characterized by rapid and straight growth with pleasant color and ease of handling [21] and has been used for a variety of purposes, such as house construction, shipbuilding, and making wooden barrels. This species could easily be vegetatively propagated, making it suitable for common garden experiments in multiple environments. Highly saturated linkage maps using a large number of genetic markers have been developed for this species [22,23], which is necessary for genetic mapping studies [13,24]. It is empirically well known that heartwood color of C. japonica could become dark or reddish in some cultivars [25][26][27]. For example, a previous study conducted in a single forest stand reported that the brightness of heartwood color (L � ) largely differs between different clones of C. japonica, and moisture and potassium content are correlated with heartwood color [28]. In addition, darker heartwood color is less preferred for the usage of non-laminated wood such as in traditional construction [28]. However, the genetic basis of genotype-by-environment interaction and the genomic regions controlling heartwood color in C. japonica is largely unknown.
This study aimed to evaluate the relative effects of genetic and environmental factors and genotype-by-environment interaction on heartwood color in C. japonica and identify genomic regions that control the heartwood color of this species using multiple common garden experiments established in three contrasting environments. These common gardens have been used to identify QTLs of wood property and growth traits and the climate sensitivity traits of C. japonica in previous studies [13,24]. This study further compared QTLs of heartwood color to previously reported QTLs of wood property, such as heartwood moisture content and density, to identify QTLs that could affect multiple heartwood property traits.

Plant material, common garden experiment, and linkage map
The plant material, common garden experiment, and linkage map previously described by Mori et al. [13,24] were used in this study. The mapping progeny comprised 139 genotypes from F1 plants, YI96 (female parent) and YI38 (male parent), were used as mapping populations in the present study. YI96 andYI38 were derived from local varieties of the Kyushu region (southern part of Japan): "Yabukuguri," "Iwao," and "Kumotooshi." Previous studies have reported differences in the heartwood color of these cultivars [25,29]. Significant differences in heartwood colors among these cultivars (P < 0.01, analysis of variance; for details, see S1 and S2 Tables) were also confirmed using wood samples of the cultivars obtained in a previous study [30]. YI96 was obtained by crossing "Yabukuguri" × "Iwao," whereas YI38 was obtained by crossing "Yabukuguri" × "Kumotooshi." Common garden experiments were established for the 139 genotypes, with three clonal replicates in three different environments in Japan: Ibaraki site (36˚11 0 4.66@N, 140˚13 0 1.80@E), Chiba site (35˚20 0 52.55@N, 1401 0 47.92@E), and Kumamoto site (32˚41 0 58.36@N, 130˚45 0 17.64@E); the experiments consisted of three clonal replicates with a randomized block design with each tree planted at an interval of 1.8 m. Local varieties 'Yanase', 'Sanbu', and 'Shakain' were planted at the surrounding edges of the study site in Ibaraki, Chiba, and Kumamoto, respectively, to minimize the edge effect [13]. Three to five soil samples from each site were collected for soil analysis to evaluate water content, total carbon (g/kg), total nitrogen (g/kg), pH (H 2 O, KCl), and exchangeable cations (Ca, Mg, K, Na; cmol(+)/kg) in 2017 (S3 Table). Two linkage maps based on 858 genetic markers (YI96: heterozygous; YI38: homozygous) and 916 genetic markers (YI38: heterozygous; YI96: homozygous) were used to construct the YI96 and YI38 linkage maps, respectively. Genetic markers with large numbers of missing values (>1%) were excluded. The average marker distance between adjacent markers in the YI96 and YI38 linkage maps was 1.47 and 1.56 cM, respectively.

Measurement of heartwood color
Wood samples from 40 to 80 cm above ground were collected from each mapping progeny to measure heartwood color. Unedged boards (5 cm) with pith were sawn using a bandsaw and dried indoors for >1 year to attain the air-dried state. The quartersawn faces of the wood samples were planed with a hand-feed planer before heartwood color measurement. Heartwood color was measured as L � a � b � color space (CIELAB; lightness L � , green to red hues a � , and blue to yellow hues b � ) and spectral reflectance (400-700 nm with 20 nm interval), which are commonly used to assess wood color in various tree species, including C. japonica [31,32]. Heartwood color traits were measured using a spectrophotometer NF333 (Nippon Denshoku Industries, Tokyo, Japan), with a measurement area and a light source of φ8 mm and D65/10˚, respectively. Heartwood color measurement was done at five equal distance points at the center and the edges (left and right) of the heartwood in a longitudinal direction, resulting in 15 measurement points per wood sample. When the heartwood width of the sample was narrow, the center of the heartwood was measured, resulting in five measurement points for those samples. For each trait and site, the observed values were averaged for each sample, and genotypes with two or more replicates per site were retained for subsequent analyses.

Phenotypic data analysis
Heartwood color can be spatially biased based on microscale environmental conditions within sites. For example, the heartwood color of C. japonica is affected by soil potassium content [33]. Thus, Moran's I test was conducted to investigate the level of spatial autocorrelation for each trait using the R package "spdep" [34]. A linear model for the block design (with block and genotype as fixed effects) was fitted for each site and trait, and the residuals were used for Moran's I tests.
Moran's I statistic ranges from −1 (perfectly dispersed) to 1 (perfectly aggregated). The significance of spatial autocorrelation was tested under the assumption of spatial randomness.
The genotypic means (best linear unbiased estimators) for each site and trait were predicted using the linear model for the block design. When there was a significant positive spatial autocorrelation with Moran's I test for a given trait within a site (P < 0.05; S4 Table), the genotypic means were predicted using the model for the block design with additional spatial autocorrelation terms (penalized spline curves; P-splines) using the R package "SpATS" [35]. The predicted genotypic means for each site were used for subsequent QTL analyses and the correlation matrix (see below).
The contribution of genetic and environmental effects to the total variability of heartwood color traits was estimated using variance components analyses. Each trait, genotype, site, block (nested within sites), and interactions between genotypes and sites were used as random variables to estimate variance components and 95% confidence intervals (CIs). Variance component analyses were performed using the "anovaVCA" function of the R package "VCA." To further understand the environmental effects of three sites on heartwood color traits, the difference in phenotypic values among the sites for each trait was tested using the post-hoc analysis of Tukey's all-pairs comparison. Phenotypic data analysis was performed using R version 4.1.1 [36], and multiple comparison tests were performed using the R package "multcomp" [37]. Correlation coefficients (Pearson's r) among the genotypic means of the heartwood color traits were obtained for each site, and the level of significance of the correlations was tested with P-values calculated by the t-distribution using the R package "Hmisc." Correlations among heartwood color traits were visualized with heatmaps of the correlation matrix for each site using the R package "corrplot" (Fig 2) and with correlation network plots using the R package "corrr" (S1 Appendix).

QTL analysis
The methodology of QTL analysis conducted in this study was fully described by Goto et al. [38]. QTL analysis was conducted using hierarchical linear models with a regularized horseshoe prior distribution. One marker from a given pair of markers that were closely located (<1 cM) and showed a high correlation (r > 0.7) with other markers was removed to reduce collinearity and redundancy among variables. The posterior distribution was estimated using four chains with 5000 iterations after the first 3000 iterations as burn-in using the R package "rstanarm" [39]. The convergence of the chains was confirmed using the Gelman-Rubin statistic (R < 1.1). The global shrinkage parameter of the regularized horseshoe prior was set using the equation proposed by Piironen and Vehtari [40]; p0 DÀ p0 1 ffi ffi n p , where p0 is the expected number of relevant variables, D is the number of variables, and n is the number of observations. p0 was set to 5 for all traits, and the sensitivity of the p0 value was confirmed by ranging p0 from 1 to 9 (S5 Table). The default values of the regularized horseshoe prior in the R package "rstanarm" were used for the remaining parameters.
The marginals of the posterior distribution of marker effects may be biased due to the high correlation among variables. Thus, a model selection method (projection-predictive variable selection) [40] was applied to the model. This method uses a model with the best predictive power (a reference model) to find a simpler model with similar predictions as the reference model (predictive projection). The linear model described above was used as the reference model. Model selection was done when variables in the reference model contained significantly more information than the null model, which was tested by leave-one-out cross-validation [41]. When the marginals of the marker effects in the selected model did not contain zero in the 95% and 80% confidence intervals (CIs), the markers were defined as significant and suggestive QTLs, respectively. The contribution of each QTL was estimated by the coefficient of determination (R 2 ) of the simple regression, and the contribution was defined as the phenotypic variation explained (PVE) for a QTL. Because the resolution of QTL positions depends on the size of the mapping progeny, QTLs that were closely located should be considered the same QTL. In this study, QTLs located within 5 cM were defined as the same QTL. QTL analysis and model selection were done using R version 4.1.1 [36]. The leave-one-out cross-validation procedure was performed using the R package "loo" [41], and the model selection method was done using the R package "projpred" [42]. Because L � a � b � color space and spectral reflectance were highly correlated (see Results), QTLs were shown only for L � , a � , and b � in the QTL analysis results to avoid complexity. The QTL analysis results of spectral reflectance can be found in the (S2 and S3 Appendices; S6-S8 Tables).

Heartwood color traits
Heartwood color traits (L � a � b � color space and spectral reflectance) of mapping progenies in the three sites are summarized in Table 1, and the sample scanned images of the sample clones (genotypes) with the highest and lowest L � are shown in Fig 1. Significant spatial autocorrelation was detected for all heartwood color traits in the Kumamoto site and a � and b � in the Ibaraki site, whereas no significant spatial autocorrelation was found for heartwood color traits in the Chiba site (S4 Table). L � in the Ibaraki site was significantly smaller (i.e., darker) than Chiba and Kumamoto sites (P < 0.05; Table 1). a � in the Chiba site was significantly larger (i.e., more reddish) than Ibaraki and Kumamoto sites (P < 0.05). This tendency was also observed in significantly higher spectral reflectance in longer wavelengths (640-700 nm) in the Chiba site (P < 0.05). b � in the Kumamoto site was significantly smaller (i.e., more blueish) than the other two sites (P < 0.05). This tendency was also observed in significantly higher spectral reflectance in shorter wavelengths (440-500 nm) in the Kumamoto site (P < 0.05).
Correlations among heartwood colors and properties (moisture content and density) are summarized in Fig 2. All heartwood colors were significantly positively correlated (P < 0.001; S9 Table); L � and spectral reflectance had a strong and significant positive correlation (r > 0.96; P < 0.001), whereas relatively small to moderate significant positive correlations were found among L � , a � , and b � (0.801 > r > 0.297; P < 0.001) at the three study sites. Heartwood moisture content was significantly negatively correlated with all heartwood color traits at the Chiba and Ibaraki sites (−0.326 > r > −0.576; P < 0.001; S9 Table) and with some heartwood color traits at the Kumamoto site (−0.180 > r > −0.220; S9 Table). Heartwood density was significantly positively correlated with a � at the three sites (0.417 > r > 0.220; P < 0.01; S9 Table).

Proportion of the variability of heartwood color traits explained by genotype and environment
The proportion of variance components of genotype, site, block, and interactions between genotypes and sites in heartwood color traits are shown in Table 2. All heartwood color traits were significantly different between genotypes and interactions between genotypes and sites in terms of 95% CI (Table 2). Genotypes, sites, and their interaction explained 14.2%-33.9%, 1.7%-14.6%, and 7.9%-15.2% of the total variability of heartwood color traits, respectively. Genotypes and environments (site and block) explained the largest and smallest variance among variables in all heartwood color traits. The variance components of genotypes were

QTL analysis
QTL analysis identified 64 significant QTLs and 12 suggestive QTLs for heartwood color traits of L � a � b � color space in the three sites (Table 3; Fig 3). The number of QTLs identified for L � , a � , and b � were 34, 21, and 23, respectively (S4 Appendix). The number of QTLs identified for

PLOS ONE
L � a � b � color space in Ibaraki, Chiba, and Kumamoto sites were 44, 13, 21, respectively. One stable QTL (a QTL identified in all sites) was identified for b � (S7 Table). Three QTLs were identified for L � a � b � color space in the two sites: one QTL was identified each in Ibaraki and Chiba sites, Ibaraki and Kumamoto sites, and Chiba and Kumamoto sites. Some of the identified QTLs were clustered in genomic regions, such as 94.0 to 99.1 cM in LG06 of the YI96 map (Fig 3). In addition, 6 of 13 QTLs with PVE > 10 were located in LG09. Two QTLs of L � were considered as the same QTLs identified for heartwood moisture content in a previous study [13] (S8 Table).

Heartwood color traits
Heartwood color traits observed for mapping progenies were comparable to previous reports on cultivars used in this study (i.e., Yabukuguri, Iwao, and Kumotooshi cultivars). For example, Kawazumi et al. [25] reported that L � of Kumotooshi and Yabukuguri cultivars were 56.8 ± 7.2 and 69.3 ± 3.0, respectively. Tsushima et al. [29] investigated the heartwood color of C. japonica cultivars in different initial spacings and reported that L � , a � , and b � were 65.2 to 69.4, 5.3 to 8.0, and 19.8 in Iwao cultivar and 72.2 to 73.7, 8.3 to 9.4, and 22.6 to 23.8 in Yabukuguri cultivar, respectively. Observations in previous reports and this study (S1 Table) on the three cultivars are mostly within the range of the observed values of L � , a � , and b � of mapping  progenies. These findings indicate that heartwood color in this species is genetically preserved in various conditions (e.g., tree age, environments), which would be important for the improvement of this traits in tree-breeding programs. Heartwood color traits, especially for L � , were negatively correlated with the heartwood moisture content, indicating that heartwood with higher moisture content is darker in color in C. japonica mapping progenies. This tendency is empirically known and has also been reported in previous studies [25,29]. For example, Tsushima et al. [29] reported that the average values of L � and green moisture content of heartwood of Yabukuguri cultivar were 72.2, 72.6, and 73.7 and 168%, 163%, and 145% in the initial spacing of 1500, 3000, and 5000 trees/ ha, respectively. In addition, Sugimoto et al. [43] reported that the lightness of C. japonica wood depends on the density and the amount of reflection in the layer, especially at long wavelengths. This study also found a weak but significant positive correlation between reflectance

PLOS ONE
at long wavelengths (>600 nm) and heartwood density at Chiba site but almost no significant correlation between reflectance at short to intermediate wavelengths (400-600 nm) and density. These results implied that moisture content is more related than density with heartwood color in C. japonica. Furthermore, these suggested that heartwood color is not only important in terms of esthetic value but also a valuable indicator that represent wood moisture which could affect drying processes of the wood [28].

Relative contribution of genotype, environment, and genotype-byenvironment interaction on heartwood color
Although the effects of genotype and genotype-by-environment interaction on heartwood color traits were all significant, the relative contribution of these variables varied largely on heartwood color traits of mapping progenies. Genetic effects alone explained 30.7% on average of the total variance of heartwood color traits, whereas environmental effects (i.e., sites) alone

PLOS ONE
only explained 3.9%, indicating that heartwood color traits are largely genetically controlled, whereas environmental differences could be minimal, which could also be suggested by the scanned images of the clones with the highest and lowest L � (Fig 1). The large contribution of genetic effects in the heartwood color of C. japonica was also reported in a previous study conducted in a single forest stand [28]. These findings indicated that heartwood colors in this species have high potential to be improved through breeding in various environmental conditions, which is particularly important for long-lived forest trees that experience substantial environmental changes such as ongoing climate change.
The effect of environmental differences on heartwood color varied among heartwood color traits. For example, environmental differences only explained 2.4% and 14.6% of the total variance in L � and a � , respectively. These results indicated that the environmental difference could affect heartwood color in terms of reddish hue (red/green) but not the brightness in C. japonica. As noted above, a � was significantly higher for mapping progenies in the Chiba site, which was the only site without significant spatial autocorrelation in either heartwood color trait. These results suggested that environmental conditions in the Chiba site could have affected the reddish hue of the heartwood of mapping progenies. Although there was not sufficient evidence to evaluate the cause of these results, the soil analysis results implied that pH and potassium content in the soil of the Chiba site was higher and lower, respectively, than the other sites. These soil compositions may have altered the reddish hue of the heartwood of mapping progenies. Indeed, previous studies reported that heartwood color in C. japonica was influenced by potassium content [33].
The interaction between genotype and environment explained 10.6% on average of the total variance and comprised 26% of the total genetic effects, indicating that genotype-by-environment interaction is important for understanding the genetic basis of heartwood color in C. japonica. Mori et al. [13] reported that the variablity of heartwood moisture content explained by genotype and genotype-by-environment interaction were 54.8% and 39.6%, respectively, using the same mapping progenies of this study, indicating that the variability explained by genotypes and genotype-by-environment interaction on heartwood moisture content was larger than heartwood color traits. These findings indicated that heartwood color traits are largely genetically controlled, whereas heartwood moisture content is more genetically controlled than heartwood color traits.

QTLs controlling heartwood color
QTL analysis identified a large number of QTLs for L � , a � , and b � , and most of these QTLs had small to moderate effects (PVE = 7.0% on average) on heartwood color traits. This tendency was consistent with previous studies that identified QTLs for wood properties and growth traits in tree species (reviewed in [1]). The number of identified QTLs differed largely between sites and traits. The largest number of significant QTLs was identified for L � , and the PVE of L � was larger than a � and b � . These indicated that brightness could be more genetically controlled by QTLs than other heartwood traits, such as a � and b � , and implied that brightness could be suitable for molecular breeding (e.g., genomic selection; [44,45]) of heartwood colors in C. japonica. In contrast, the largest and smallest numbers of QTLs were identified in Ibaraki and Chiba sites, respectively. As noted above, the soil conditions of the Chiba site could be substantially different compared to other sites, suggesting that the environmental conditions in the Chiba site may have hindered the expression of QTLs.
One QTL identified for b � was detected in the three sites. In addition to this QTL, three QTLs identified for spectral reflectance of short wavelengths (<520 nm) were also detected in the three sites (S7 Table). These QTLs were considered stable to environmental changes. In contrast, most of the identified QTLs were detected uniquely in either site, indicating that QTLs are mostly affected by environments, consistent with a previous study of the mapping progeny [13]. The stable QTLs identified for b � were located within the QTL cluster of 94.0 to 99.1 cM in LG06 of the YI96 map. In addition, seven of eight QTLs identified for L � , a � , and b � located within this QTL cluster were detected in two or more sites, and these QTLs had relatively larger PVE (9.3% on average) than the other QTLs (6.7% on average). These findings suggested that there might be some important genes that control heartwood color traits in a variety of environments near the QTL cluster.
Some QTLs had pleiotropic effects on heartwood color and moisture content. For example, a QTL (m_1540) located in LG09 of the YI38 map was detected in heartwood color traits of L � and spectral reflectance of 400-700 nm (S6 Table) and this QTL was also identified for heartwood moisture content in the same study site ( [13]; S6 Table). This might indicate that there could be some important pleiotropic genes controlling both heartwood moisture content and color in C. japonica near the QTL. This further indicates the importance of identifying QTLs not only for a single trait of interest but also for multiple traits in understanding the genetic basis of the complex traits in forest trees.

Conclusion
This study evaluated the relative contribution of genetic and environmental effects using multiple common gardens in contrasting environments and identified QTLs controlling heartwood color traits in C. japonica. Overall, genetic effects, including genotype-by-environment interaction, had a large effect on heartwood color traits (L � a � b � color space and spectral reflectance), whereas environmental effect alone explained less variance. Therefore, not only genetic effects alone but also genetic variation in phenotypic plasticity were considered important for breeding of heartwood color in C. japonica under changing environments. These findings showed that heartwood color in C. japonica has high potential and is suitable in terms of improving complex traits through breeding programs in wide range of environmental conditions. In addition, this study found a large number of QTLs with relatively small effects; however, some QTLs were expressed stably under multiple environments or had pleiotropic effects on heartwood color and moisture content. These findings indicated that molecular breeding approaches such as MAS and genomic selection are possible in improvement programs on heartwood color in C. japonica; MAS may be possible in the improvement programs conducted within a small range of environmental conditions, while further studies on identification of genes controlling heartwood color are essential for successful marker-based improvements. On the other hand, genomic selection would be effective in improvement programs under wide range of environmental conditions. These necessitates further studies with different genomic approaches such as association mapping studies [4,46,47] to validate the findings of this study and to capture wider genetic variations because the present study was limited to segregated populations derived from local cultivars.
Supporting information S1